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The moment of core collapse in star clusters with a mass 
function 
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ABSTRACT 

Star clusters with multi-mass components dynamically evolve faster than those mod- 
elled with equal-mass components. Using a series of direct A^-body simulations, we in- 
vestigate the dynamical evolution of star clusters with mass functions, especially their 
core collapse time. Multi-mass clusters tend to behave like systems with a smaller 
number of particles, which we call the effective number of particles (iVoff) and for 
which NcS = M/rriraax (here M and rrimax are the total cluster mass and the mass 
of the most massive star in the cluster, respectively). We find that the time of core 
collapse is inversely proportional to the mass of the most massive star in the cluster 
and analytically confirm that this is because the core collapse of clusters with a mass 
function proceeds on the dynamical friction timescale of the most massive stars. As 
the mass of the most massive star increases, however, the core-collapse time, which is 
observed as a core bounce of the cluster core from the evolution of the core density or 
core radius, becomes ambiguous. We find that in that case the total binding energy 
of the hard binaries gives a good diagnosis for determining the moment at which the 
cluster core collapses. Based on the results of our simulations, we argue that the core 
bounce becomes ambiguous when the mass of the most massive star exceeds 0.1% of 
the total mass of the cluster. 

Key words: galaxies: star clusters: general — methods: numerical 



1 INTRODUCTION 

Star clusters are collisional systems with a negative heat 
capacity, and therefore they dynamically evolve to eventu- 
ally reach core collapse. The process to core collapse is sim- 
ply described using a semi-analytic treatment of the energy 
tran sfer from the cluster core to the outer part of the clus- 
ter IjLvnden-Bell fc Woodlll968l ). and it was confirmed us- 
ing various methods such as g aseous models (JNakada 1978; 
iBettwieser &: Sugimotd Il984l). Fokker-Planck simulation s 



I lnagakilll980l : ICohnlll980l : iLvnden-Bell fc Eggletonlll98d) 
and direct A/"-body simu lations l|Bettwieser &: Sugimotd 
1 19851 : iMakino et all Il986l '). Because of the high stellar 
density of the core durin g its collapse , dynamical bina- 
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ries form in this phase (jSpitzer fc Hard Il97ll : IA, 
1 1974 ). Once binaries form in the core, they gener- 
ate energy which is transporte d outward by interaction s 
with other cluster members (|Heggid 1 19751 : iHutI 1 19831 ). 



Due to the binary heating the core finally bounces, 
and then the core o scillates when the number o f stars 
is sufficiently large (|Sugimoto fc Bettwieseij 119831 ). Such 
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grav othermal o scillation was first found using gaseous mod- 
els! Suginioto_& Bcttwicscr 1983; Bcttwicscr & Sugimotd 
1 1984 ) but the behaviour was later confir med using direct A*'- 
bodv JMakino et al.llT98^ : lMakindl 19961 ) and Fokker-Planck 
(|Goodmanlll987l ) simulations. 

The time between cluster birth and the moment of core 
collapse (what we will call the core-collapse time, or tec) 
is proportional to the two-bo dy relaxation time at the half- 
mass radius trh (|Spitzeilll987l ). For clusters in which all stars 
have the same mass, the core collapses at tec — 15irh. This 
result has a theoretical background and was confirmed with 
simulatio ns using a wide variety of techniques: Monte-Carlo 
methods llSpitzer fc Shullll 19751 ). Fokker-Plan ck calculation s 
l|Cohdll979 ^. and direct JV-bodv simulations (|Makindll996l ). 
With a spectrum of stellar masses, however, the core- 
collapse time be comes much shorter than w hen all stars have 
the same mass Ijlnagaki fc Wivantdll984 ). From direct A'^- 
body simulations, it is obtained that tec — 0.2trh with a 
realistic mass function IjPortegies Zwart fc McMillanll2002l ). 

ICiirkan et al.l l|2004 ) performed a large number of 
Monte Carlo simulations covering a wide range of the con- 
centration parameters and mass functions. They showed 
that the core-collapse time scales with the central relaxation 
time, ire rather than with irh. Here the central relaxation 
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time is defined as 



"'"~ G2(m)pc In A' ^^-^ 

where (Jc,3d and pc are the central three-dime nsional veloc- 
ity d ispersion and core density, respectively (|Spitzerlll987l : 
[Gurkan et al. 2004). The Coulomb logarithm, In A, is writ- 
ten as ln7A'^, where A'^ is the total number of stars. For 
the central relaxation tim e of star clusters, it is n umerically 
obtained that 7 = 0.015 (JGiersz fc HeggieHlQQd ) (the clas- 
sic theoretical value for t^h with single-mass components is 
-Y = 0.11 Spitzen (il987i )). 

In the case with a power-law mass function in which 
m-max and (m) are the maximu m mass and th e mean mass 
of the stellar mass distribution, iGiirkan et al.l (2004) found 
that the core-collapse time tec oc {mmux/ (m))^^"^ . They 
also demonstrated the existence of a minimum to the core- 
co llapse time, whi c h is t cc/^rc = 0.15. This value is also seen 
in lGoswami et al.l (|2012i ). in whic h they use the s ame M onte 
Carlo code that was adopted bv lCiirkan et al.l (|200J), but 
with a wider range of initial conditions including initially 
mass-segregated models. The arguments for the particular 
exponent (being —1.3) and the minimum to the core-collapse 
time, however, was not discussed and remains unclear. 

In this paper we show the results of core collapse simu- 
lations of star clusters with power-law mass functions using 
direct A'^-body simulations. We find that the core collapse 
time scal es tec oc (mm a. x/i'm) )'^ , contrary to the earlier 
finding of Giirkan et al.l l|2004h . but we support our finding 
with analytic arguments. We further argue that core col- 
lapse is driven by the sinking of the most massive star to 
th e cluster centre, by dynamica l frictio n (as was suggested 
in iPortegies Zwart fc McMillanI (|2002l ')'). The minimum to 
the core-collapse time then corresponds to the time required 
for the most massive star to reach the cluster centre. For 
the most extreme mass functions the minimum core-collapse 
time then naturally depends on the crossing time of the sys- 
tem. The core bounce becomes less pronounced for larger 
values of m-mn^/ (jn), and this is qualitatively understood 
from the dynamical evolution being driven by the most mas- 
sive stars in the cluster. For a mass ratio Af/77imax ~ 10'^ core 
collapse becomes hard to determine, and it even becomes in- 
distinguishable for Af/rrimax ~ 100, because in those cases 
the core will eventually be composed of only a few massive 
stars, almost irrespective of the total number of stars in the 
cluster. 



2 A^-BODY SIMULATIONS AND THE INITIAL 
CONDITIONS 

We perf ormed a se ries of A^-body simulations using King 
models (|Kinel 119661 ) with a non-dimensional concentration 
parameter. Wo, of 3 and 6 as initial density profiles. Here- 
after, we adopt A'-body units in which, G = M = —AE — 1, 
where G, M, and E are the gravitational constant, the to- 
tal mass, and the total energy of the cluster, respectively 
l|Heggiel ll975l lH In table [1] we summarize the initial con- 
ditions for the runs. For each simulation we adopted A' 
particles from N = 2048 (2k), 8192 (8k), 32768 (32k) to 

^ see also |http : //en ■ wlklpedla ■ org/wikl/M-body_miits| 



A^ = 1310172 (128k), with a power-law mass function with 
exponent a and an upper- mass limit of mmax- The value 
of ?nniax/(?Ti) is Varied from 1.0 (equal mass) to 517, but 
for models with fewer particles we adopted a smaller value 
of nimax- For the mass function exponent a, we adopted 
a = 2.35 (|SalpeteJll955h . 1.7, and 1.2. In Tables [2] and E] we 
summarize the simulation results. 

All simulations are performed us ing the sixth-order 
predi ctor-corrector Hermite scheme (JNitadori fc Making 
l2008n running on GPU using the Sapporo2 lib rary 
(Bell eman et al.ll2008l : iBedorf fc Portegies Zwartll2012l ) and 
also on CPU cl usters using the two-d imensional paralleliza- 
tion sche me by Nitadori et al.l (120061 ). We used a time step 
criterion (jNitadori fc Making l2008h with accuracy parame- 
ter, 77 =0.1-0.3. The energy error ^ 10~^ for all simulations 
over the entire duration of the simulation. In some cases the 
time step of the calculations dropped below ~ 10~^^, which 
occurs near the moment of the formation of a binary of a few 
lOOfcT for the models and with muiax/{m) ^ 2. (We express 
binding energies in terms of kT = {m)afY>, where ctid is the 
one-dimensional velocity dispersion of the cluster, 3NkT/2 
is the initial kinetic energy of the entire system.) To prevent 
even smaller time steps, which have unpleasant consequences 
for the performance, we adopted a small softening e; for the 
simulations with Wo = 3 we adopted e = 1/(200 A') and 
e = 1/(130A') for models with Wo = 6. With this softening 
we are able to resolve binaries with a semi-major axis of a 
1200fcT and 780fcT for the simulations with Wo = 3 and 
Wo = 6, respectively. 



3 CORE COLLAPSE IN MULTI-MASS 
CLUSTERS 

The core-collapse time is determined by the moment of core 
bounce, which is seen in the time evolution of the core radius 
or core density. For some models it is difficult to distinguish 
the core collapse, because there does not seem to be a peak 
in the density evolution or a depression in the core radius. 
Another indicator for determining the moment of core col- 
lapse is by monitoring the evolution of the binding energy of 
dynamically formed binaries. During the core collapse, hard 
binaries form in the cluster. They are hardened by three- 
body encounters in the cluster core and generate the energy 
for the core bounce. In this section, we present the evolu- 
tion of the core radius and the density from the simulations, 
and then we discuss the relation between the core evolution 
and the dynamically formed binaries in order to provide an 
objective determination of the moment of core collapse and 
ultimately define the core collapse. 



3.1 The evolution of core density and radius 

In Figure [T] we present the evolution of the core density 
and the core radius for an equal-mass model with Wo ~ 
3, A^ =2k (left) and the same model but with a Salpeter 
mass function with mma.yi/{m) = 8 (right) . We calculate the 
core r adius and density using a method of lCasertano fc Hull 
(|l985r ), but we took into account the mass of the central 
particles. Compared to the equal mass case, the core collapse 
in the models with mmax/(?7i) ;^ 8 is less clear, although the 
core noticeably expands after a slight depression. 
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Table 1. Properties of King Models 



Wo 


Pc 


CTclD 


CTlD 


re 


fh 


Mc 


3 
6 


0.652 
2.11 


0.518 
0.503 


0.520 
0.636 


0.543 
0.293 


0.839 
0.804 


0.238 
0.117 



If we increase m^^^/ {m), the core collapse becomes 
more ambiguous. In Figure [2] we present the results of the 
same simulation as we presented in Figure [l] but with 
mmax/{m) = 32 (left in Fig.© and 129 (right). In these 
simulations the mass of the most massive star is comparable 
to the mass of the initial cluster core. After the core shrinks, 
it keeps the small core radius and slowly expands (see right 
panel of Figure [2]) . 

Based on our simulations we conclude that a core 
bounce occurs if M/mmax ~ 10'^, and the rapid core expan- 
sion becomes apparent only if M/mmax ~ 100. We illustrate 
these in the left panel of Figure[3l where we present the same 
models as in Figure [2] but with A'' = 128k. This model satis- 
fies the first criterion (M/mmax ~ 10'^) and as expected the 
core bounce is clearly visible. With 7Tiniax/{w,) = 517 (the 
right panel of Figure [S]), however, which does not satisfy ei- 
ther of the criteria and as a consequence the core bounce 
is obs cured. These criteria are similar to the Spitzer insta- 
bility l|ST3itzeilll987l ). In the case of multi-mass components, 
the criterion for the Sp itzer instability is M/rrimax ~ 10'' 
l|Breen fc Heggiell2012al ). We discuss these criteria further 
in section 4. 



3.2 The determination of core-collapse time using 
hard binaries 

In the previous subsection we discussed the lack of a core 



bounce for the case where M/m„ 



10 . In particular if 



Af/rrimax ^ 100, the core in these cases however still ex- 
pands quite dramatically after some time. In these cases it 
becomes very hard to use the core size or density peak to 
determine the moment of core collapse, but the expansion of 
the core indicates that something like a core collapse must 
have happened. In order to quantify this moment we fo- 
cus on the evolution of hard (dynamically formed) binaries, 
which are suspected to generate the energy that causes the 
core to bounce. 

In the bottom panels in Figures [U [21 and [31 we present 
the total binding energy of hard binaries in the various simu- 
lations. We accumulated the binding energy of hard binaries 
(with > 3kT), but the total binding energy is almost com- 
parable to that of the hardest binary. By comparison of the 
evolution of the binding energy with the core density (or core 
radius), we observe that the moment of the core bounce is 
consistently occurring at the same moment that the total 
binding energy in hard binaries reaches a few hundreds of 
kT. 

We estimate the critical binding energy for the core 
bounce from a discussion on the energy emitted by the hard 
binary via a three-body encounter. In the dynamical evo- 
lution of star clusters through core collapse, the cluster re- 
sponds to a core collapse by a bounce, and the occurrence is 
associated with the moment when the energy produced by 



har d binaries e xceeds the potential energy of the cluster core 
(po (|Hutlll996l : fHeggie fc Hutll2003l '). The energy released in 
an encounter between a single star and a binary with bind- 
ing energy E^in is ASbin = 0.4£;bin (|Heggielll975l '). The core 
then bounces when the binding energy of the hard binary 
satisfies 

iSbin > o^("^><^0. (2) 

Following the discussion in iHutI (|l996l l , assuming that the 
core is virialized at the moment of the bounce, the potential 
energy of the core is (f>o = 2K = Nc{'m)a^ 3J3. Here we adopt 
o"c,3D for the velocity dispersion of the stars in the core and 
K is the kinetic energy of the core. In order to be able to 
bounce the core, the binaries should have a binding energy 
of at least 



Ebin 



— A"c(m)a,%D 



(3) 

(4) 



~ 200{m)<3D. 

Here we adopted A^c ~ 80 l|Hutlll996l ). Although ctcsd 
is the velocity dispersion at the moment of core bounce, 
we assume that it is similar to the initial core veloc- 
ity dispersion, (Tc,o,3d, because the core evolves foUow- 
ing pc PC r"" and k = 2.2-2.3 during the core collapse 
fCoimlgS^; 'Lvnde n-Bell fc Eggletoiilll980l : lTakahashi|[l995l : 
Heggie fc Hut 2003i ), and the refore the velocity dis persion 
changes only by a factor of 2 (JGiersz fc Heggielll996l '). From 
the definition of kT, the critical binding energy for core 
bounce is obtained as i?bin,cr ~ (200cr^^o,3i3/'^iD)fc7^ ~ 
{600a^Qi]j/aiYy)kT. In the case of King models with Wo = 3 
and 6, iJbin.cr ~ QOOkT and 390fcT, respectively. This es- 
timate is consistent with the peak in the binding energy 
distribution of escaping binaries obtained from a direct A^- 
b ody simulations in equal-mas s simulations (see Figure 22 
in lTanikawa fc Fukushigell2009l ) . In Figure [4] we present the 
relation between the core-collapse time obtained from the 
evolution of core density (or radius) and that obtained from 
the earlier discussed requirement Ehin > Ehin,cT- Both ap- 
pear to predict consistently the same moment for the core 
collapse, and we can adopt the binary binding energy argu- 
ment to obtain an objective and physically motivated argu- 
ment for the moment of core bounces, and therefore for the 
moment of core collapse. Hereafter we adopt this criterion 
to define the moment of core collapse. 



3.3 The core-collapse time 

In Figure [SI we present tcc/irc as a function of mmax/(m). 
Here we defined the core-collapse time tec as the moment 
when the total biding energy of hard binaries in the clus- 
ter exceeds QOOkT and 390fcT for models with Wo ^ 3 
and Wo = 6, respectively. The core-collapse time for single- 
mass component models is tcc/trc — 20 for Wo ~ 3 and 
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Table 2. Models with Wo = 3 



Model 


Wo 


a 


mmin/(ni) 


»nniax/(m> 


N 


e 


tec 


teookT 


w3-2k-eq 


3 


- 


1.0 


1.0 


2k 


l/(200Af) 


735 


750 


w3-8k-eq 


3 


- 


1.0 


1.0 


8k 


l/(200Af) 


2.41 X 10^ 


2.41 X 10^ 


w3-2k-m2-Sal 


3 


2.35 


0.607 


2.02 


2k 


l/(200Af) 


435 


442 


w3-8k-m2-Sal 


3 


2.35 


0.607 


2.02 


8k 


l/(200Af) 


1.26 X 10^ 


1.26 X 10^ 


w3-2k-m8-Sal 


3 


2.35 


0.391 


8.07 


2k 





85 


92 


w3-8k-m8-Sal 


3 


2.35 


0.391 


8.07 


8k 





297 


300 


w3-2k-m32-Sal 


3 


2.35 


0.329 


32.3 


2k 





33 


39.5 


w3-8k-m32-Sal 


3 


2.35 


0.329 


32.3 


8k 





55 


61.2 


w3-32k-m32-Sal 


3 


2.35 


0.329 


32.3 


32k 





217 


225 


w3-2k-ml29-Sal-l 


3 


2.35 


0.296 


129.2 


2k 





13 


12.6 


w3-2k-ml29-Sal-2 


3 


2.35 


0.296 


129.2 


2k 





19 


22.1 


w3-8k-ml29-Sal 


3 


2.35 


0.296 


129.2 


8k 





35 


29.4 


w3-32k-ml29-Sal 


3 


2.35 


0.296 


129.2 


32k 





74.6 


81.9 


w3-128k-ml29-Sal 


3 


2.35 


0.279 


516.6 


128k 





175 


172 


w3-2k-m258-Sal-l 


3 


2.35 


0.286 


258.3 


2k 





_ 


15.6 


w3-2k-m258-Sal-2 


3 


2.35 


0.286 


258.3 


2k 





- 


19.3 


w3-8k-m258-Sal 


3 


2.35 


0.286 


258.3 


8k 





22.8 


21.1 


w3-32k-m258-Sal 


3 


2.35 


0.286 


258.3 


32k 





62 


64.0 


w3-2k-m517-Sal-l 


3 


2.35 


0.279 


516.6 


2k 





_ 


15.6 


w3-2k-m517-Sal-2 


3 


2.35 


0.279 


516.6 


2k 





- 


22.9 


w3-8k-m517-Sal 


3 


2.35 


0.279 


516.6 


8k 





- 


23.9 


w3-32k-m517-Sal 


3 


2.35 


0.279 


516.6 


32k 





30.0 


30.1 


w3-128k-m517-Sal 


3 


2.35 


0.279 


516.6 


128k 





40.7 


36.3 


w3-2k-m2-al.2 


3 


1.2 


0.442 


2.02 


2k 


l/(200Ar) 


341 


346 


w3-8k-m2-al.2 


3 


1.2 


0.442 


2.02 


8k 


l/(200Af) 


1.19 X 10^ 


1.14 X 10^ 


w3-2k-m8-al.2 


3 


1.2 


0.0336 


8.07 


2k 





88 


127 


w3-8k-m8-al.2 


3 


1.2 


0.0336 


8.07 


8k 





294 


294 


w3-2k-m32-al.2 


3 


1.2 


5.28 X 10-4 


32.3 


2k 





32 


39.3 


w3-8k-m32-al.2 


3 


1.2 


5.28 X 10-4 


32.3 


8k 





99 


96.8 


w3-2k-ml29-al.2 


3 


1.2 


5.16 X 10-'' 


129.2 


2k 





_ 


18 


w3-8k-ml29-al.2 


3 


1.2 


5.16 X 10-6 


129.2 


8k 





- 


26.5 


w3-32k-ml29-al.2 


3 


1.2 


5.16 X lO-f' 


129.2 


32k 





92.8 


98.6 


w3-2k-m258-al.2 


3 


1.2 


2.13 X 10-* 


258.3 


2k 





_ 


17.8 


w3-8k-m258-al.2 


3 


1.2 


2.13 X 10-* 


258.3 


8k 





- 


16.4 


w3-8k-m517-al.2 


3 


1.2 


1.23 X 10-8 


516.6 


8k 





9.7 


8.94 


w3-32k-m517-al.2 


3 


1.2 


1.23 X 10-* 


516.6 


32k 





- 


46.1 


w3-2k-m2-al.7 


3 


1.7 


0.525 


2.02 


2k 


l/(200Ar) 


414 


426 


w3-2k-m8-al.7 


3 


1.7 


0.192 


8.07 


2k 





98 


108 


w3-2k-m32-al.7 


3 


1.7 


0.0856 


32.3 


2k 





36 


37.6 


w3-8k-m32-al.7 


3 


1.7 


0.0856 


32.3 


8k 





70.5 


73.3 


w3-2k-ml29-al.7 


3 


1.7 


0.0403 


129.2 


2k 





11 


10.8 


w3-8k-m258-al.7 


3 


1.7 


0.0304 


258.3 


8k 





- 


11.6 



tec is the core-collapse time obtained from the core density evolution, while tsooi^-p is the time in which the total binding energy of hard 
binaries in the cluster exceed 600fcT. Here we defined teookT (for models with Wo = 3) and tsgokT (for models with Wo = 6) as the 

time when, for the first time, Bbin > ^bin cr- 
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Wq = 3, A' = 2k, equal mass 




M/q = 3, N = 2k, nj„^Tx/(m) = 8 




Figure 1. Time evolution of the cluster parameters. To panel: the evolution of the core density. Middle panel: evolution of the core 
radius (solid) and the 10%, 30, 50 and 80% (bottom to top) Lagrangian radii (dotted curves). Bottom panel: the total binding energy in 
hard binaries. The left three panels give the results for the model with Wq = 3 for equal mass stars, and the right three panels give the 
results for a mass function with a = 2.35 and minax/(("i) = 8. 



10-* 

103 

100 

10 

1 



Wa = 3, jV = 2k, ??J„^Tx/(nj) = 32 












, iiriuJAjUy 


vw 


JM^jJi-t^^ 







Wo = 3,N = 2k, m,,,J{m) = 129 





Figure 2. Same as Figure[T] but for models with mma^K / {{in) = 32 (left) and 129 (right). 
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Table 3. Models with Wo = 6 



Model 


Wo 


a 


mmiii/(m) 


m.max/(»n) 


N 


e 


tec 


i390kT 


w6-2k-eq 


6 


- 


1.0 


1.0 


2k 


l/(130Af) 


650 


652 


w6-8k-eq 


6 


- 


1.0 


1.0 


8k 


l/(130Af) 


1.86x10^ 


1.86x10^ 


w6-2k-m2-Sal 


6 


2.35 


0.607 


2.02 


2k 


l/{130Af) 


230 


251 


w6-8k-m2-Sal 


6 


2.35 


0.607 


2.02 


8k 


l/{130Af) 


806 


808 


w6-2k-m8-Sal 


6 


2.35 


0.391 


8.07 


2k 





51.8 


52.9 


w6-8k-m8-Sal 


6 


2.35 


0.391 


8.07 


8k 





114 


138 


w6-2k-m32-Sal 


6 


2.35 


0.329 


32.3 


2k 





19.9 


13.1 


w6-8k-m32-Sal 


6 


2.35 


0.329 


32.3 


8k 





31.9 


31.9 


w6-2k-ml29-Sal 


6 


2.35 


0.296 


129.2 


2k 





15.5 


17.3 


w6-8k-ml29-Sal 


6 


2.35 


0.296 


129.2 


8k 





18.6 


15.4 


w6-32k-ml29-Sal 


6 


2.35 


0.296 


129.2 


32k 





34 


35.1 


w6-2k-m517-Sal 
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50 for Wq = 6, whi ch are consistent with previous results 
(JGiirkan et al.ll2004l . and references therein). With a mass 
function, tcc/trc decreases as we increase myaax/{m), and it 
follows (mmax/("i))~^ (dashed line in Figure [Sj as far as 
mmax/(m) < 30. 

We analytically derive the core-collapse time for mod- 
els with stellar mass functions using the dynamical friction 
timescale of the most massive stars in the cluster, assum- 
ing that star clusters collapse when the most massive stars 
reach the cluster centre. The dynamical friction timescale of 
the most massive star with mmax is esti mated from a simple 
equat ion. We follow the description in lBinnev fc Tremaind 
l|200a ) in section 8.1.1., in which the dynamical friction 
timescale of a black hole which spirals in to the galactic cen- 
tre is derived. We assume that the massive star has initially 
a circular orbit with velocity, Vc, at a distance r from the 
cluster centre. From equation (8.9) in iBinnev fc Tremaind 
1(2003), the frictional force F = mui3.^\dvm / dt\ on the mas- 
sive star is 



47rG2mL,p(r)lnA' 



erf(X) 



2X 



exp(— X 



(5) 



dimensional velocity dispersion and Vc is equivalent to 
the two-dimensional velocity dispersion. The value of the 
Coulomb logarithm here is different from t hat in equation 
( [ill, an d therefore we describe it as In A' (see lSpinnato et al.l 
(j2003r )). In order to simplify this equation, we assume the 
density distribution to be a singular isothermal sphere, 
pir) — v'^ / {i-rrGr^) , and then equation ((S| becomes 



0.428 In A' 



/ G^^max 



(6) 



The angular momentum change of the massive star due to 
the friction is 



dL 
dt 



= ^Fr 



0.428 In A' 



,Gm 



2 

max 



(7) 
(8) 



In an isothermal sphere the circular velocity is independent 
of radius, and the angular momentum at radius r is written 
as L = m^^^rVc. With equation Q, we obtain 



r—- = —0.428 . 

dt Vc 



(9) 



where X 



c/(\/2(73d) ~ 1, where ctsd is the three- The dynamical friction timescale of the most massive star is 
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finally written as 
1.65 r^CTsD 



tdt 



(10) 



In A' Gmirn^ ' 

This result is shown in Figure [S] (dashed curves) and it is 
consistent with the simulations. H ere, we adopted r as the 
half-mass radius and A' = O.IA^ (|Giersz fc Heggid [l994al : 



iHeggJe fc Hu3l2003l ). For the Coulomb logarithm for Uc, 
we adopted A' = O.lA'^c, where A^'c = McN is the number 
of particles in the core. The dynamical friction timescale 
of the most massive star is proportional to m^l_^, and as 
a consequence the core-relaxation time is proportional to 
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(m) ^ (see Eq.[T}, and therefore we obtain that fcc/irc oc 

(mniax/(m))"-^. 



4 DISCUSSION 



3.4 The minimum core-collapse time 

As we see in Figure O models with fewer particles (A'' — 2k) 
deviate from the analytic result when rriynax/ (rn) ~ 30. The 
critical value of nimax/ (m) , however, depends on N. For 
models with A'^ = 128k the simulations and theory give 
consistent results until m-mn^/ (rn) ~ 500. The models with 
N between 2k and 128k show a consistent picture of the 
moment that the models deviate from the theory for an 
increasing value of 7Timax/('Ti) with respect to A^. We bol- 
ster our earlier conclusion that a core bounce requires that 
Af/mmax ^ 100. The core-collapse time saturates at smaller 
mmax/im) for models with fewer particles. We will discuss 
this critical value of M/rrimax in section [TTI and we will now 
focus to estimating the "minimum" core-collapse time we 
observe in Figure O 

The minimum to the core-collapse time depends on A', 
which is a result of the dependency of Uc on A^. We con- 
sider that this minimum core-collapse time depends on the 
crossing time, icross, of the cluster, because the dynamical 
friction time cannot be shorter than the crossing time. We 



adopt tc 



)^"/(7iD, where (r^)^" is the root mean 



square of the distance from the cluster centre. The minimum 
core-collapse time obtained from the simulations is roughly 
consistent with lOf cross- The dotted lines in Figure [S] give 
lOtcross. Those lines also depend on N, because f cross does 
not depend on A', whereas ire does. 



4.1 



A-dependence and comparison with 
single-component models 



!t 10^ and for rapid core expansion in 
^ 100. Here we will make an analogy 



We discuss the criterion for an objective core collapse in the 
cases where M/rrim ^ ^ '^'^ 
the case that M/m 
with single-component models. From a wide range of ana- 
lytic calculations and simulations, it is well established that 
the dynamical evolution of star clusters such as relaxation, 
core collapse, core bounce, and gravothermal oscillation de- 
pends on the total number of particles in the system (and in 
the core). For example, gravothermal oscill ation occurs only 
when the num ber of particles exceed ~ 10* (|Goodmanlll987l : 
iMakindl 19961 ). This criterion comes from the number of par- 
ticles in the core after core bounce, Ac- (Here we define Ac 
as the averaged number of particles after core collapse, and 
we adopt Acb as the number of particles at the moment of 
core bounce.) The gravothermal oscillation occurs only when 
Ac > Acb. While the va lue of Acb is considered a constant 
in th e range of 10 to 80 (|Goodmanlll987l : lMakind[T99^: IHuI 
19961), A'c depends on A^. If we adopt Ac ~ A^^^ JMakinc 



1996). we can confirm that A > 10* satisfies Ac > Acb. 

The behaviour after core collapse changes when N de- 
creases. So long as A'c ~ A^cb (i.e., 10^ ;S A" < 10*) we 
observe similar evolution but the collapse becomes shal- 
lower for smaller A^, and the gravothermal oscillations dis- 
appear (see Fig. 1 in [Makino 1996). For 10 < A < 10^ the 
core b ounce disappears (see Figure 10 in iGiersz fc Heggiel 
ll994b|). These transitions are quite similar to those we ob- 
serve if we increase mmax- Our results of the multi-mass case 
are scalable to those of the single mass case if we define an 
effective number of particles (Aoff) in the latter. For an equal 
mass system, Aoff = A', but when we introduce a spectrum 
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Figure 6. The evolution of the density profile for model w3-2k- 
ml29-Sal-2 (the same model in the right panel of Figure [2]l. 



of masses, A'eff = M/m,max. Interestingly, a similar conclu- 
sion is obtained from rece nt results for two -component and 
multi-component systems IjBreen fc Heggiell2012bl lari. 

If the number of particles drops below ~ 100, the be- 
haviour of the A^-body system changes from being a many- 
body system to a few-body system, which evolves chaotically 
IjMcMillan et al.l [l988) rather than deterministic. In few- 
body systems it is hard to notice any core collapse in the 
evolution of core density and radius. In these cases the bi- 
nary cannot harden to ^ lOOfcT because the total energy 
budget of the cluster ^ lOOkT. In such systems the hard 
binary stops interacting with other cluster members when it 
becomes too tight and the surrounding density becomes too 
low. As a re sult, the "t enured" or "frozen binary" remains in 
the cluster (|Casertand ll985). This situation is seen in multi- 
mass model if Af/mmax ~ 100 (see Figure[2]). The dynamical 
evolution of frozen binaries almost stops in this case, but the 
frozen binary is still in the cluster, and likely to be in the 
core. In Figure [6]we present the evolution of density profiles 
of model w3-2k-ml29-Sal-2 (we used the same model in the 
right panel of Figure [2]). The frozen binary and its relatively 
low-density environment are noticeable as the high density 
peak in the centre and a depletion in the more extended 
core. This effect is similar to the core mass-deficiency ar- 
guments used in the nuclei with bin ary black holes after a 
major galaxy merger (|Merrittll2010l ). In the formation pro- 
cess of frozen binaries, massive particles concentrated in the 
cluster centre are ejected from the cluster by sling-shot in- 
teractions with the binary. This mechanism leads to the for- 
mation of massive runaway stars arou nd dense, young star 
clusters (JFuiii fc Portegies Zwartll20lll '). 



4.2 Dynamical evolution driven by the most 
massive stars 

The dynamical evolution of equal-mass models scales with 
trc irrespective of the particle number (|Spitzerlll987l ). We 
confirm this hypothesis with our multi-mass simulations, so 
long as Af/nimax /t 100. For multi-mass models, however, 
icc/irc decreases for increasing m^s.^/ {rn){^ /max), and we 
demonstrate that ice is determined by the dynamical friction 
timescale of the most massive stars in the cluster. Here we 
show that the dynamical evolution of star clusters with a 
mass function is driven by the most massive stars in the 
cluster. 

The relaxation time of a multi-mass system is dom- 
inated by the dynamics of the most massive star, 
and the global relaxation time is a factor of F-m = 
ln(7Af//niax)/(/max ln(7A'^)) shortcr than that of an equal 
mass system. In Figure [T] we present the evolution of the 
core density scaled in time by the product of frc and F^i. 
As long as the model satisfies the condition of core collapse 



{M/mn 



10 ), the evolutionary tracks of the core den- 



sity are scaled with FmUc- The scaling parameters are deter- 
mined using the initial cluster. When the cluster approaches 
core collapse, the models with a mass function start to de- 
viate from the equal-mass case. 

The maximum core density, pc,max, reached the moment 
of the core collapse (the depth of the core collapse) decreases 
as mina.x/{m} increases. We present the relation between 
Pc,max and muiEux/im) in Figured This phenomenon is sim- 
ilar to the relation betw een the maximum density and A'' 
IjGiersz fc Heggiell994bl). If we t ry to understand this follow- 
ing the discussion in iHud l| 19961 ) , in which it is derived that 
p oc A''-^, we expect that the maximum density decreases 
following (mmax/(wT))-^ assuming that the core consists of 
the most massive stars and therefore that p oc N~^. In Fig- 
ure [8l however, the power is shallower than —2, although 
we see the trend that the maximum density decreases for 
increasing m^lax/{n^). This might be because the mean par- 
ticle mass in the core is smaller than rrimax. Note that we 
can measure the core density only in a snapshot, which is 
limited by our output frequency. We therefore are likely to 
miss the highest density peak; it is very difficult to catch the 
exact moment of the highest density in an A'^-body simula- 
tions. This can be solved by storing the particle position and 
velocity inf ormation in a time resolved data format , as w as 
proposed in lFaber et"all l|2010l ) (see also lFarr et~aLl l|2012l )'). 



4.3 Comparison with previous results 

We find that tcc/irc oc (?nmax /(?7t))~^, which is inconsis- 
tent with a previous work by iGiirkan et al.l (|2004), who 
conclude that the core-collapse time is proportional to 
(?nmax/("^))-^'^. This discrepancy can be attributed to their 
use of a Monte-Carlo code, which cannot properly follow the 
cluster all the way to the moment of core collapse, but shows 
core bounce at some earlier instance. A hint of the early ter- 
mination of their calculations is visible in their Fig. 2, where 
the Lagrangian radii of the simulations with a mass func- 
tion suddenly truncates. They identify this moment as core 
collapse, but by inspection of our own simulations the sud- 
den break in the inner most Lagrangian radii is generally 
associated with the formation of the first hard (few kT) 
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binary. We tested this hypothesis by analysing the resuhs 
of our simulation up to the moment of the formation of 
the first hard {3kT) binary. In figure [9] we present the for- 
mation time of the first 3fcT-binary, t^kT, for models with 
Wo = 3 as a function of mmax/ (rn). With our adopted for- 
mation of the first hard binary our results become consis- 



Figure 9. The moment of the formati on of the first 3kT bi- 
nary (the hard binary limit according to iTanikawa fc Fukushigd 
l2009t) as a function of the maximum mass of the cluster parti- 
cles, mmax/(m) for Wo = 3 and Salpeter MF models. Models 
with minax/('7i) = 517 initially host 3kT binaries, and there- 
fore they are not plotted in this figure. Dashed line shows 20 X 
('Timax/(ni))~"''^. The dotted lines indicate the takT/^rc = 0.15, 
which is suggested the minimum ratio between i^c and i^c by 
iGurkan et"allll2004) . 



tent with those of iGiirkan et al.l (|2004l ) . and we argue that 
t-ikr/trc oc (m,nax/{m))~-''^ and saturates at tcc/irc = 0.15. 



5 CONCLUSIONS 

We performed a series of A'^-body simulations of star clus- 
ters with various mass ranges and power laws of the mass 
function, and found that the core-collapse time follows 

tcc/trc OC {minax/{m))~^ for cluSterS with M/Tllmax ^ 100. 

When M/rrimax ~ 100, this relation breaks and ice saturates 
at ~ lOicross . Using analytic arguments we derived that star 
clusters with a mass function reach core collapse in the dy- 
namical friction timescale of the most massive stars. We also 
showed that the dynamical evolution of star clusters with a 
mass function are driven by the relaxation timescale of the 
most massive stars. We define an effective number of stars 
A^eff = Af/mmax for which a multi-mass cluster shows a sim- 
ilar core collapse behaviour as in the equal-mass case. 

When the mass of the most massive stars is relatively 



small {M/m-n 



10 ), we notice a pronounced peak in the 



core density during the evolution. As we increase the mass of 
the most massive stars, rrimax, the peak at the core bounce 
becomes lower. We found that the binding energy of the 
hard binaries are a good indicator for detecting the moment 
of the core collapse, even if the density peak is ambiguous. 
We adopted that the critical binding energy of hard bina- 
ries, with which the binary emits sufficient energy to bounce 
the core and eject stars from the system, as the moment of 
the core bounce. We conclude that the binding energy crite- 
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rion gives a more robust indicator for the core-collapse time 
compared to inspection of the core density and radius. 
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